#include "nr.h"
/* usando o teorema da compensacao */

void calcula_jacobiano1(int nelements,
                        int nelectrodes,
                        int nnodes,
                        int CURR_PATTERN,
                        double **volt_pattern_matrix,
                        double corr,
                        int no_terra,
                        double **k_global,
                        double ***k_locais,
                        int *vec_electrodes,
                        double *distancias_centroides,
                        int nnp_local,
                        double *v_rho_0,
                        double ***H,
                        double *volt_calculado_diff,
                        int *VARIABLE_TABLE,
                        int *erro_saida)
{
  int i,m,n;
  int q_mais,q_menos;
  double *current;
  double *aux_line;
  double **k_global_inversa;
  double **Hj_diff;
  double **Hj_se;
  double **volt_calculado_SE,*volt_calculado_temp;
  
  
  volt_calculado_SE=dmatrix(1,nnodes,1,2); /* a segunda coluna nao eh utilizada*/

  k_global_inversa=dmatrix(1,nnodes,1,nnodes);
  aux_line=dvector(1,nnodes);
  
  for (m=1;m<=nnodes;m++){
    aux_line[m] = k_global[m][no_terra];
    k_global[no_terra][m]=0.0;
    k_global[m][no_terra]=0.0;
  }
  k_global[no_terra][no_terra]=1.0;
  
  
  for (i=1;i<=nelectrodes;i++){
    
    volt_calculado_SE=dmatrix(1,nnodes,1,2); /* a segunda coluna nao eh utilizada*/
    Hj_se=dmatrix(1,nelectrodes,1,nelements);
    
    current=dvector(1,nnodes);
    q_mais=i;
    q_menos=(q_mais + CURR_PATTERN)%nelectrodes+1;
    /*      q_terra=(q_mais + 1 + CURR_PATTERN)%nelectrodes+1;*/
    current[vec_electrodes[q_mais]]= corr;  
    current[vec_electrodes[q_menos]]= -corr;
  
  /* resolve o problema direto */
    
    sist_lin_esparso_intel(nnodes,
                           k_global,
                           1,
                           current,
                           volt_calculado_SE,
                           erro_saida,
                           1);
    
    for (j=1;j<=nelectrodes;j++){
      for (k=1;k<=nelectrodes;k++){
      
      }
    }

    free_dmatrix(Hj_se,1,nelectrodes,1,nelements);
    free_dmatrix(volt_calculado_SE,1,nnodes,1,2); /* a segunda coluna nao eh utilizada*/
    free_dvector(current,1,nnodes);
  }
      

  invert_intel(nnodes,
	       k_global,
	       k_global_inversa,
	       erro_saida,
	       1);  /* 1: matriz simetrica   2:matriz nao simetrica*/

  /**********      restaura a matriz global          ********************/
  for (m=1;m<=nnodes;m++)
    k_global[no_terra][m] = k_global[m][no_terra] = aux_line[m];

  if((*erro_saida)!=0){
    free_dmatrix(k_global_inversa,1,nnodes,1,nnodes);
    free_dvector(aux_line,1,nnodes);
    (*erro_saida)=501;
    return;
  }

  volt_calculado_SE=dvector(1,nelectrodes);
  volt_calculado_temp=dvector(1,nelectrodes);
  Hj_diff=dmatrix(1,nelectrodes,1,nelements);
  (*H)=dmatrix(1,nelectrodes*nelectrodes,1,nelements);              VARIABLE_TABLE[20]=1;
    
  for (i=1;i<=nelectrodes;i++){
    Hj_se=dmatrix(1,nelectrodes,1,nelements);
    current=dvector(1,nnodes);

    q_mais=i;
    q_menos=(q_mais + CURR_PATTERN)%nelectrodes+1;
    /*      q_terra=(q_mais + 1 + CURR_PATTERN)%nelectrodes+1;*/
    current[vec_electrodes[q_mais]]= corr;  
    current[vec_electrodes[q_menos]]= -corr;

    monta_sensibilidade (Hj_se,
			 k_global_inversa,
			 k_locais,
			 vec_electrodes,
			 current,
			 volt_calculado_SE, /* voltagens do problema direto SE*/
			 nnodes,
			 nelements,
			 nelectrodes,
			 nnp_local,
			 distancias_centroides, /**/
			 v_rho_0);

    /*******************   Calcula Hj_diff=volt_pattern_matrix*Hj_se ********************/

    SEtoDIF(volt_pattern_matrix,
	    Hj_se,
	    nelectrodes,
	    nelements,
	    Hj_diff);

    /*******************   calcula potencial calculado differencial ********************/

    multiplyV(nelectrodes,
	      nelectrodes,
	      &(volt_pattern_matrix[1][1]),  /*  Entrar com &(Matrix1[1][1]) ou &(Matrix1[0][0]) */
	      &(volt_calculado_SE[1]),  /*  Entrar com &(Vector1[1]) ou &(Vector1[0]) */
              NULL,
	      &(volt_calculado_temp[1]),  /*  Entrar com &(Vector2[1]) ou &(Vector2[0]) */
	      22);

    for (m=1;m<=nelectrodes;m++){
      for (n=1;n<=nelements;n++)
	(*H)[((nelectrodes)*(i-1))+m][n]=Hj_diff[m][n];
      volt_calculado_diff[nelectrodes*(i-1)+m]=volt_calculado_temp[m];
    }

    free_dmatrix(Hj_se,1,nelectrodes,1,nelements);
    free_dvector(current,1,nnodes);
  }
      
  free_dvector(volt_calculado_SE,1,nelectrodes);
  free_dvector(volt_calculado_temp,1,nelectrodes);
  free_dmatrix(k_global_inversa,1,nnodes,1,nnodes);
  free_dvector(aux_line,1,nnodes);
  free_dmatrix(Hj_diff,1,nelectrodes,1,nelements);
}
